rm(list=ls())
library(metafor)

###################
#### Loading Data and AK Estimators
###################
PATH="C:/Users/hong_/Desktop/World Development/FinalCodes_Data/"
DT<-as.data.frame(read.csv(paste(PATH, "data_sorted.csv", sep='')))
source(paste(PATH, 'AK_Estimators.R', sep=''))


###################
#### Creating Dataset and attache to R working directory
###################
RegData <- as.data.frame(cbind(effect=DT$pcc, se=DT$sepcc, DT))
attach(RegData)


###################
#### Making Table 4 - Upper Panel (Estimation Results)
###################
Table4 <- as.data.frame(matrix(0, ncol=8, nrow=12))
colnames(Table4) <- c('FE', 'FEwSE', 'RE', 'REwSE', 'AK1', 'AK1wSE', 'AK2', 'AK2wSE')
rownames(Table4) <- c('mu', '95CI', 'Obs', 'LLK', 'Parameters', 'AIC', 'BIC', 'tau', 'Test1', 'p(Test1)', 'Test2', 'p(Test2)')


## Column 1 (FE)
reg_FE <- summary(rma.uni(pcc~1 ,vi=(sepcc^2),method='FE'))
reg_FE <- robust(reg_FE, id)
Table4[,1] <- as.data.frame(t(cbind(reg_FE$beta[1], paste('(', round(reg_FE$ci.lb,4), ',', round(reg_FE$ci.ub,4),')', sep=''),
                                    reg_FE$k.all, reg_FE$fit.stats[1,1],
                                    1, reg_FE$fit.stats[3,1], reg_FE$fit.stats[4,1], 0, 0,0,0,0)))
## Column 2 (FEwSE)
reg_FEwSE <- summary(rma.uni(pcc~ 1 + sepcc, vi=(sepcc^2),method='FE'))
reg_FEwSE <- robust(reg_FEwSE, id)
Table4[,2] <- as.data.frame(t(cbind(reg_FEwSE$beta[1], paste('(', round(reg_FEwSE$ci.lb[1],4), ',', round(reg_FEwSE$ci.ub[1],4),')', sep=''),
                                    reg_FEwSE$k.all, reg_FEwSE$fit.stats[1,1],
                                    2, reg_FEwSE$fit.stats[3,1], reg_FEwSE$fit.stats[4,1], 0, 0,0,0,0)))
## Column 3 (RE)
reg_RE <- summary(rma.uni(pcc~1 ,vi=(sepcc^2),method='ML'))
reg_RE <- robust(reg_RE, id)
Table4[,3] <- as.data.frame(t(cbind(reg_RE$beta[1], paste('(', round(reg_RE$ci.lb,4), ',', round(reg_RE$ci.ub,4),')', sep=''),
                                    reg_RE$k.all, reg_RE$fit.stats[1,1],
                                    2, reg_RE$fit.stats[3,1], reg_RE$fit.stats[4,1], sqrt(reg_RE$tau2), 0,0,0,0)))
## Column 4 (REwSE)
reg_REwSE <- summary(rma.uni(pcc~ 1 + sepcc, vi=(sepcc^2),method='ML'))
reg_REwSE <- robust(reg_REwSE, id)
Table4[,4] <- as.data.frame(t(cbind(reg_REwSE$beta[1], paste('(', round(reg_REwSE$ci.lb[1],4), ',', round(reg_REwSE$ci.ub[1],4),')', sep=''),
                                    reg_REwSE$k.all, reg_REwSE$fit.stats[1,1],
                                    3, reg_REwSE$fit.stats[3,1], reg_REwSE$fit.stats[4,1], sqrt(reg_REwSE$tau2), 0,0,0,0)))
## Column 5 (AK1)
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1))
reg_AK1 <- AK1Estimator(AKdata)
ci <- reg_AK1$EstResults[1,3]+qt(c(0.025,0.975), (reg_AK1$SampleSize-ncol(reg_AK1$EstResults)))*reg_AK1$EstResults[2,3]
Table4[,5] <- as.data.frame(t(cbind(reg_AK1$EstResults[1,3], paste('(', round(ci[1],4), ',', round(ci[2],4),')', sep=''),
                                    reg_AK1$SampleSize, reg_AK1$LogLiklihood,
                                    ncol(reg_AK1$EstResults), reg_AK1$aic, reg_AK1$bic, reg_AK1$EstResults[1,1], 0,0,0,0 )))
## Column 6 (AK1wSE)
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1, se=sepcc))
reg_AK1wSE <- AK1Estimator(AKdata)
ci <- reg_AK1wSE$EstResults[1,3]+qt(c(0.025,0.975), (reg_AK1wSE$SampleSize-ncol(reg_AK1wSE$EstResults)))*reg_AK1wSE$EstResults[2,3]
Table4[,6] <- as.data.frame(t(cbind(reg_AK1wSE$EstResults[1,3], paste('(', round(ci[1],4), ',', round(ci[2],4),')', sep=''),
                                    reg_AK1wSE$SampleSize, reg_AK1wSE$LogLiklihood,
                                    ncol(reg_AK1wSE$EstResults), reg_AK1wSE$aic, reg_AK1wSE$bic, reg_AK1wSE$EstResults[1,1], 0,0,0,0 )))
## Column 7 (AK2)
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1))
reg_AK2 <- AK2Estimator(AKdata)
ci <- reg_AK2$EstResults[1,5]+qt(c(0.025,0.975), (reg_AK2$SampleSize-ncol(reg_AK2$EstResults)))*reg_AK2$EstResults[2,5]
Table4[,7] <- as.data.frame(t(cbind(reg_AK2$EstResults[1,5], paste('(', round(ci[1],4), ',', round(ci[2],4),')', sep=''),
                                    reg_AK2$SampleSize, reg_AK2$LogLiklihood,
                                    ncol(reg_AK2$EstResults), reg_AK2$aic, reg_AK2$bic, reg_AK2$EstResults[1,1], 0,0,0,0 )))
## Column 8 (AK2wSE)
AKdata <- as.data.frame(cbind(id, effect=pcc, se=sepcc, constant=1, se=sepcc))
reg_AK2wSE <- AK2Estimator(AKdata)
ci <- reg_AK2wSE$EstResults[1,5]+qt(c(0.025,0.975), (reg_AK2wSE$SampleSize-ncol(reg_AK2wSE$EstResults)))*reg_AK2wSE$EstResults[2,5]
Table4[,8] <- as.data.frame(t(cbind(reg_AK2wSE$EstResults[1,5], paste('(', round(ci[1],4), ',', round(ci[2],4),')', sep=''),
                                    reg_AK2wSE$SampleSize, reg_AK2wSE$LogLiklihood,
                                    ncol(reg_AK2wSE$EstResults), reg_AK2wSE$aic, reg_AK2wSE$bic, reg_AK2wSE$EstResults[1,1], 0,0,0,0 )))
Table4 <- as.matrix(Table4)


###################
#### Making Table 4 - Lower Panel (Test Results)
###################
## Column 1 (FE)
TestPair <- c(2,1) # Column 2 (unrestricted) vs Column 1 (Restricted)
test1.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
TestPair <- c(3,1)
test2.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
Table4[9:12,1] <- as.matrix(rbind('FE(w/SE)', test1.p.value,'RE', test2.p.value))

## Column 2 (FEwSE)
TestPair <- c(4,2) 
test1.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
Table4[9:12,2] <- as.matrix(rbind('RE(w/SE)', test1.p.value,'NA', 'NA'))

## Column 3 (RE)
TestPair <- c(4,3)
test1.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
TestPair <- c(5,3)
test2.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
Table4[9:12,3] <- as.matrix(rbind('RE(w/SE)', test1.p.value,'AK1', test2.p.value))

## Column 4 (REwSE)
TestPair <- c(6,4)
test1.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
Table4[9:12,4] <- as.matrix(rbind('AK1wSE', test1.p.value,'NA', 'NA'))

## Column 5 (AK1)
TestPair <- c(6,5)
test1.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
TestPair <- c(7,5)
test2.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
Table4[9:12,5] <- as.matrix(rbind('RE(w/SE)', test1.p.value,'AK2', test2.p.value))

## Column 6 (AK1wSE)
TestPair <- c(8,6)
test1.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
Table4[9:12,6] <- as.matrix(rbind('AK2wSE', test1.p.value,'NA', 'NA'))

## Column 7 (AK2)
TestPair <- c(8,7)
test1.p.value <- pchisq(2*(as.numeric(as.matrix(Table4[4,TestPair[1]]))-as.numeric(as.matrix(Table4[4,TestPair[2]]))), (as.numeric(as.matrix(Table4[5,TestPair[1]]))-as.numeric(as.matrix(Table4[5,TestPair[2]]))), lower.tail=FALSE)
Table4[9:12,7] <- as.matrix(rbind('AK2wSE', test1.p.value,'NA', 'NA'))

## Column 8 (AK2wSE)
Table4[9:12,8] <- as.matrix(rbind('NA', 'NA','NA', 'NA'))

Table4 <- as.data.frame(Table4)
View(Table4)
detach(RegData)
write.csv(as.data.frame(Table4), "EstimationResults/Table04.csv", row.names=TRUE)